home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / cos.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  6.9 KB  |  311 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    cos   double precision cosine
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    cos
  29.  *    machine independent routines
  30.  *    trigonometric functions
  31.  *    math libraries
  32.  *
  33.  *  DESCRIPTION
  34.  *
  35.  *    Returns double precision cosine of double precision
  36.  *    floating point argument.
  37.  *
  38.  *  USAGE
  39.  *
  40.  *    double cos (x)
  41.  *    double x;
  42.  *
  43.  *  REFERENCES
  44.  *
  45.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  46.  *    1968, pp. 112-120.
  47.  *
  48.  *  RESTRICTIONS
  49.  *
  50.  *    The sin and cos routines are interactive in the sense that
  51.  *    in the process of reducing the argument to the range -PI/4
  52.  *    to PI/4, each may call the other.  Ultimately one or the
  53.  *    other uses a polynomial approximation on the reduced
  54.  *    argument.  The sin approximation has a maximum relative error
  55.  *    of 10**(-17.59) and the cos approximation has a maximum
  56.  *    relative error of 10**(-16.18).
  57.  *
  58.  *    These error bounds assume exact arithmetic
  59.  *    in the polynomial evaluation.  Additional rounding and
  60.  *    truncation errors may occur as the argument is reduced
  61.  *    to the range over which the polynomial approximation
  62.  *    is valid, and as the polynomial is evaluated using
  63.  *    finite-precision arithmetic.
  64.  *    
  65.  *  PROGRAMMER
  66.  *
  67.  *    Fred Fish
  68.  *
  69.  *    68881 support included by Michael Ritzert
  70.  *
  71.  *  INTERNALS
  72.  *
  73.  *    Computes cos(x) from:
  74.  *
  75.  *        (1)    Reduce argument x to range -PI to PI.
  76.  *
  77.  *        (2)    If x > PI/2 then call cos recursively
  78.  *            using relation cos(x) = -cos(x - PI).
  79.  *
  80.  *        (3)    If x < -PI/2 then call cos recursively
  81.  *            using relation cos(x) = -cos(x + PI).
  82.  *
  83.  *        (4)    If x > PI/4 then call sin using
  84.  *            relation cos(x) = sin(PI/2 - x).
  85.  *
  86.  *        (5)    If x < -PI/4 then call cos using
  87.  *            relation cos(x) = sin(PI/2 + x).
  88.  *
  89.  *        (6)    If x would cause underflow in approx
  90.  *            evaluation arithmetic then return
  91.  *            sqrt(1.0 - x**2).
  92.  *
  93.  *        (7)    By now x has been reduced to range
  94.  *            -PI/4 to PI/4 and the approximation
  95.  *            from HART pg. 119 can be used:
  96.  *
  97.  *            cos(x) = ( p(y) / q(y) )
  98.  *            Where:
  99.  *
  100.  *            y    =    x * (4/PI)
  101.  *
  102.  *            p(y) =  SUM [ Pj * (y**(2*j)) ]
  103.  *            over j = {0,1,2,3}
  104.  *
  105.  *            q(y) =  SUM [ Qj * (y**(2*j)) ]
  106.  *            over j = {0,1,2,3}
  107.  *
  108.  *            P0   =    0.12905394659037374438571854e+7
  109.  *            P1   =    -0.3745670391572320471032359e+6
  110.  *            P2   =    0.134323009865390842853673e+5
  111.  *            P3   =    -0.112314508233409330923e+3
  112.  *            Q0   =    0.12905394659037373590295914e+7
  113.  *            Q1   =    0.234677731072458350524124e+5
  114.  *            Q2   =    0.2096951819672630628621e+3
  115.  *            Q3   =    1.0000...
  116.  *            (coefficients from HART table #3843 pg 244)
  117.  *
  118.  *
  119.  *    **** NOTE ****    The range reduction relations used in
  120.  *    this routine depend on the final approximation being valid
  121.  *    over the negative argument range in addition to the positive
  122.  *    argument range.  The particular approximation chosen from
  123.  *    HART satisfies this requirement, although not explicitly
  124.  *    stated in the text.  This may not be true of other
  125.  *    approximations given in the reference.
  126.  *            
  127.  */
  128.  
  129. # include <stdio.h>
  130. # include <math.h>
  131. # include "pml.h"
  132.  
  133. #if !defined (__M68881__) && !defined (sfp004)    /* mjr++        */
  134.  
  135. static char funcname[] = "cos";
  136.  
  137. static double cos_pcoeffs[] = {
  138.     0.12905394659037374438e7,
  139.    -0.37456703915723204710e6,
  140.     0.13432300986539084285e5,
  141.    -0.11231450823340933092e3
  142. };
  143.  
  144. static double cos_qcoeffs[] = {
  145.     0.12905394659037373590e7,
  146.     0.23467773107245835052e5,
  147.     0.20969518196726306286e3,
  148.     1.0
  149. };
  150.  
  151.  
  152.  
  153.     double y;
  154.     double yt2;
  155.     struct exception xcpt;
  156.  
  157. double cos (x)
  158. double x;
  159. {
  160.     if (x < -PI || x > PI) {
  161.     x = fmod (x, TWOPI);
  162.         if (x > PI) {
  163.         x = x - TWOPI;
  164.         } else if (x < -PI) {
  165.         x = x + TWOPI;
  166.         }
  167.     }
  168.     if (x > HALFPI) {
  169.     xcpt.retval = -(cos (x - PI));
  170.     } else if (x < -HALFPI) {
  171.     xcpt.retval = -(cos (x + PI));
  172.     } else if (x > FOURTHPI) {
  173.     xcpt.retval = sin (HALFPI - x);
  174.     } else if (x < -FOURTHPI) {
  175.     xcpt.retval = sin (HALFPI + x);
  176.     } else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) {
  177.     xcpt.retval = sqrt (1.0 - (x * x));
  178.     } else {
  179.     y = x / FOURTHPI;
  180.     yt2 = y * y;
  181.     xcpt.retval = poly (3, cos_pcoeffs, yt2) / poly (3, cos_qcoeffs, yt2);
  182.     }
  183.     return (xcpt.retval);
  184. }
  185.  
  186. #endif    /* !__M68881 && !sfp004    */
  187. #ifdef    sfp004
  188.  
  189. __asm("
  190.  
  191. comm =     -6
  192. resp =    -16
  193. zahl =      0
  194.  
  195. ");    /* end asm    */
  196.  
  197. #endif    sfp004
  198. #if defined (__M68881__) || defined (sfp004)
  199.     __asm(".text; .even");
  200.  
  201. # ifdef    ERROR_CHECK
  202.  
  203.     __asm("
  204.  
  205. _Overflow:
  206.     .ascii \"OVERFLOW\\0\"
  207. _Domain:
  208.     .ascii \"DOMAIN\\0\"
  209. _Error_String:
  210.     .ascii \"cos: %s error\\n\\0\"
  211. .even
  212. | pml compatible cosgent
  213. | m.ritzert 7.12.1991
  214. | ritzert@dfg.dbp.de
  215. |
  216. |    /* NAN  = {7fffffff,ffffffff}        */
  217. |    /* +Inf = {7ff00000,00000000}        */
  218. |    /* -Inf = {fff00000,00000000}        */
  219. |    /* MAX_D= {7fee42d1,30773b76}        */
  220. |    /* MIN_D= {ffee42d1,30773b76}        */
  221.  
  222. .even
  223. double_max:
  224.     .long    0x7fee42d1
  225.     .long    0x30273b76
  226. double_min:
  227.     .long    0xffee42d1
  228.     .long    0x30273b76
  229. NaN:
  230.     .long    0x7fffffff
  231.     .long    0xffffffff
  232. p_Inf:
  233.     .long    0x7ff00000
  234.     .long    0x00000000
  235. m_Inf:
  236.     .long    0xfff00000
  237.     .long    0x00000000
  238. ");
  239. #endif    ERROR_CHECK
  240.  
  241. __asm("
  242. .even
  243.     .globl _cos
  244. _cos:
  245. ");    /* end asm    */
  246.  
  247. #endif    /* __M68881__ || sfp004    */
  248. #ifdef    __M68881__
  249.  
  250.     __asm("
  251.     fcosd    a7@(4), fp0    | cos
  252.     fmoved    fp0,a7@-    | push result
  253.     moveml    a7@+,d0-d1    | return_value
  254.     ");    /* end asm    */
  255.  
  256. #endif    __M68881__
  257. #ifdef    sfp004
  258.     __asm("
  259.     lea    0xfffa50,a0
  260.     movew    #0x541d,a0@(comm)    | specify function
  261.     cmpiw    #0x8900,a0@(resp)    | check
  262.     movel    a7@(4),a0@        | load arg_hi
  263.     movel    a7@(8),a0@        | load arg_low
  264.     movew    #0x7400,a0@(comm)    | result to d0
  265.     .long    0x0c688900, 0xfff067f8    | wait
  266.     movel    a0@,d0
  267.     movel    a0@,d1
  268.     ");    /* end asm    */
  269.  
  270. #endif    sfp004
  271. #if defined (__M68881__) || defined (sfp004)
  272. # ifdef    ERROR_CHECK
  273.     __asm("
  274.     lea    double_max,a0    |
  275.     swap    d0        | exponent into lower word
  276.     cmpw    a0@(16),d0    | == NaN ?
  277.     beq    error_nan    |
  278.     swap    d0        | result ok,
  279.     rts            | restore d0
  280. ");
  281. #ifndef    __MSHORT__
  282. __asm("
  283. error_nan:
  284.     moveml    a0@(24),d0-d1    | result = +inf
  285.     moveml    d0-d1,a7@-
  286.     movel    #62,_errno    | NAN => errno = EDOM
  287.     pea    _Domain        | for printf
  288. ");
  289. #else    __MSHORT__
  290. __asm("
  291. error_nan:
  292.     moveml    a0@(24),d0-d1    | result = +inf
  293.     moveml    d0-d1,a7@-
  294.     movew    #62,_errno    | NAN => errno = EDOM
  295.     pea    _Domain        | for printf
  296. ");
  297. #endif    __MSHORT__
  298. __asm("
  299. error_exit:
  300.     pea    _Error_String    |
  301.     pea    __iob+52    |
  302.     jbsr    _fprintf    |
  303.     addl    #12,a7        |
  304.     moveml    a7@+,d0-d1
  305.     rts
  306.     ");
  307. # else    ERROR_CHECK
  308. __asm("rts");
  309. # endif    ERROR_CHECK
  310. #endif /* __M68881__ || sfp004    */
  311.